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We revisit the model of a Brownian particle in a heat bath submitted to an actively controlled 
force proportional to the velocity that leads to thermal noise reduction {cold damping). We inves- 
tigate the influence of the continuous feedback on the fluctuations of the total entropy production 
and show that the explicit expression of the detailed fluctuation theorem involves different dynamics 
and observables in the forward and backward processes. As an illustration, we study the analyti- 
cally solvable case of a harmonic oscillator and calculate the characteristic function of the entropy 
production in a nonequilibrium steady state. We then determine the corresponding large deviation 
function which results from an unusual interplay between 'boundary' and 'bulk' contributions. 
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I. INTRODUCTION 



In parallel with the recent developments in nanotechnology and single-molecule manipulations, there is an increasing 
interest in understanding the stochastic energetics of small systems driven away from thermal equilibrium. In this 
context, fluctuation theorems (FT) play a central role as they describe exact symmetry properties of the probability 
distributions of various thermodynamic quantities such as work, heat, or entropy (see Refs. [IHlj for recent reviews). 
Typically, a (detailed) FT is a relation of the form 



PjSt = S) 
P{St = -5*) 



(1) 



where St is an observable integrated over a time interval t and C is a positive constant. In other words, a FT states 
that positive fluctuations are exponentially more probable than negative ones, which can be generally ascribed to the 
breaking of time-reversal symmetry at the level of stochastic trajectories. 

Recently, the extension of stochastic thermodynamics and fluctuation theorems to systems under feedback control 
has become an active field of research (see Ref. [5] and references therein). Feedback loops, in which some microscopic 
information about the state of the system is used to manipulate its evolution, are indeed important in many engineering 
applications and also play a crucial role in biological motors. Feedback control may be interpreted as a kind of 
"Maxwell's demon" |6^, which requires to generalize the second law of thermodynamics and to modify the various 
fluctuation theorems. 

One of the simplest example of a small classical system under a continuous feedback control is a Brownian "particle" 
in a heat bath submitted to a velocity-dependent external force. This results in a reduction of thermal fluctuations, 
as illustrated for instance by experiments with an atomic force microscope ( AFM) [71 [8] . This technique is named 
"cold damping" and is now used in a wide variety of optomechanical or electromechanical systems (see e.g. Ref. [5] 
for a review on optomechanical cooling and Refs. [101 E] for an application to a gravitational wave detector). A 
theoretical description of such a (classical) molecular refrigerator was provided in Refs. [12 US] where it was shown 
that the contraction of phase space induced by the additional viscous damping force could be interpreted as an entropy 
pumping mechanism. In particular, it was claimed in Ref. [T3] that the FT takes the form of Eq. (with ( = 1) 
when the entropy pumping term is included in the overall entropy production. The goal of the present study is 
to reexamine this statement within the path-integral formalism of Langevin dynamics and to specify the dynamics 
and observables that are associated to the probabilities appearing in the numerator and denominator of Eq. ([I]). 
There is indeed a subtlety in the definition of the so-called "backward" (or time-reversed) trajectory in a system 
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with velocity-dependent feedback, which makes the observables measured during the forward and backward processes 
different. Ahhough the issue of time-reversal symmetry of the feedback force was careful discussed in Ref.[12|, this 
point was not clearly presented in Ref. [13] and not fully appreciated in the subsequent literature. Here, for the sake 
of simplicity, we mainly consider the case of a time-independent nonequilibrium steady state (NESS). We also assume 
that the measurement process is error free, in contrast with most recent works on feedback control which investigate 
the properties of the mutual information acquired through a discrete series of measurements [5l I14j . Our discussion is 
illustrated by the paradigmatic (but experimentally relevant) example of a harmonic oscillator for which calculations 
can be performed analytically. 

The paper is organized as follows. In section 2, we briefly recall the analysis performed in Refs. [12 [13] and 
we (re)-derive the integral and detailed fluctuation theorems for the entropy production in the specific case of a 
feedback control proportional to the velocity, which corresponds to the actual experimental situation [3 |H1 [IHl fTT] . 
The stochastic harmonic oscillator is studied in section 3 where we exactly determine the characteristic function of 
the entropy production in a steady state and investigate the properties of the large deviation function (the complete 
asymptotic form of the probability distribution function is given in appendix). We conclude in section 4. 



II. ENTROPY PRODUCTION AND FLUCTUATION THEOREMS 



We consider a single Brownian particle (or "system" ) in contact with a heat bath at temperature T whose dynamics 
is governed by the one-dimensional underdamped Langevin eauation|13| 

mvs = -jvs - d^^Va{xs) + g{vs) + (2) 

where Vg = is is the velocity of the particle at time s, 7 is the friction coefficient and is a delta-correlated white 
noise with variance 2jT (Boltzmann's constant is set to unity throughout this work). Va{x) is a potential that can be 
externally controlled via a time-dependent parameter a{t), and g{v) is a velocity-dependent force that results from a 
feedback mechanism which detects the motion of the particle in real time, like in the AFM experiments described in 
Refs. [HIH]. The force g{v) is a source of entropy production and at constant a the system eventually reaches a NESS 
where heat is permanently dissipated. The probability distribution p{x, v, t) of the system in phase space is solution 
of the corresponding Kramers equation 

dtp{x, V, t) = -vdxp{x, V, t) -f —d^ I [^v + dxVa{x) - g{v)]p{x, v, t) + j—dyp{x, v,t)\ . (3) 

TO TO J 

Let us first briefly recall the analysis of Refs. [T^Hn] for the entropy production along a single trajectory {xs}se[o,t] 
of the system during a time interval < s < t. Within the stochastic energetics (or thermodynamics) framework [31 [2D], 
one usually identifies two contributions to the entropy production: 

i) the entropy change in the medium, which corresponds to the heat dissipated in the environment, 

1 1 /■* 

ASm[{xs}] = 7^Q[{xs}] = - I ds Xsi-yxs - 6] 







1 

T 



ds Xs[mxs + dx^Va{xs) - givs)] , (4) 



where the sign of Q[{a:s}] is here chosen to be positive if the heat flows out of the system into the heat bath, 
ii) the entropy change in the system itself[2Tl [22] 



AS = - lnp{xt,vt,at) + \iip{xo,VQ,ao) , (5) 

where the probability distribution, solution of Eq. is evaluated along the stochastic trajectory. (Note that all 
products of stochastic quantities as in Eq. Q are defined with the Stratonovich prescription [20].) Using the Kramers 
equation, one can show that a third contribution appears in the presence of the velocity-dependent force g{v), 

1 /■* 

ASpu[{xs}] = — / dsdy^g{vs) , (6) 
m Jo 

which is interpreted in Refs. |12lll3j as an "entropy pumping" performed by the external agent that manipulates the 
feedback force (this may be for instance an optical or electromechanical device). In the case of a friction-like control. 
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this contribution is negative (see below Eq. 
only the combination 



M). All these quantities fluctuate from one trajectory to another and 



(7) 



J] = AS^[{xs}] +AS- ASpu[{xs}] 



is always non- negative when performing the ensemble average jl2j. £[{2;^}] can thus be interpreted as the overall 
entropy production in the "super-system" composed of the particle, the heat bath and the external agent. On 
the other hand, the ensemble average of AStot = ASm + AS, the entropy production in the particle (or system) 
and the bath, can be negative. The velocity-dependent feedback thus implies a modification of the second law of 
thermodynamics. 

In order to discuss the fluctuation theorems in a more specific framework, we now assume that g{v) is proportional 
(but opposite) to the particle velocity, i.e. g{v) = —j'v (with 7' > 0), like in the AFM setup described in 0[S]. 
Thermal fluctuations are reduced by this additional friction force and the effective temperature of the system in 
a steady state, defined by its average kinetic energy, becomes lower than the heat bath temperature (hence heat 
permanently flows from the bath to the system). Specifically, in the case of a harmonic potential, one has T^ff = 
771 < >— r7/(7 + 7')[71|51[TD] (see also section 3 below). For a linear feedback, the entropy pumping contribution 
defined by Eq. ^ does not depend on the stochastic trajectory and it decreases linearly with the observation time 



ASp-a — 



-t 



(8) 



It is well known that one can relate AS'm[{a;s}], the entropy change in the medium, to the ratio of the probability 
functional for the forward and backward (i.e. time-reversed) trajectories[lT]. Taking into account the presence of 
the additional friction coefficient 7', the probability of the path {xs}s£[o,t]j given that the system started in the state 
{xo,vo), has the following expression 



V[{xs}\xo,vq] = Cexp 



7 + 7' 
2m 



1 



477 



ds (mXs -I- (7 + 7')is + dj:Va{Xs)J 



(9) 



where C is a normalization factor (see Ref. [23 for an explicit derivation of the path probability associated to an 
underdamped Langevin equation with a general non-conservative force). In the present case, it is crucial not to 
include the linear term (7 -t- ^')t/ {2m) appearing in the exponential into the normalization factor. Indeed, since the 
probability of the path {x{s)} s£[o t] defined by the time-reversal operation Xs = Xt-s, Vg = —vt-s, Ois = at-s, is given 
by 



V[{xs}\xo,vo] = C exp 



1 + 1' 
2m 



t- 



1 



47r 



ds (^mxs - (7 + l')^s + dxVaixs)^ 



(10) 



where {xo,vo) = {xt,—vt), one must also change the sign of 7' in order to extract AS'm[x(s)] from the ratio of the 
two probabilities. This yields 



'P+[{xs}\xo,vo] 

V-[{Xs}\xo,Vo] 



exp 



/ dsXs(mXi,+-f'xs + d^Va{xs)) 
m T Jo \ ) 



exp{AS'm[{a;s}] - ASpu} 



(11) 



where the subscripts -I- and — refer to the trajectories obtained with 7' and —7', respectively. Choosing the appropriate 
backward path associated with a given forward path is always an issue in a nonequilibrium state. One has indeed the 
choice between changing or not changing the sign of the external parameters that specify the state, and the proper 
choice is the one that leads to a physically meaningful result for the concrete system under considcration[24 . This is 
the case here, but one must emphasize that changing 7' into —7' is not a benign transformation: thermal fluctuations 
are then enhanced instead of being damped, and the Langevin dynamics does not lead to a stationary state at constant 
a if the effective friction coefficient 7 — 7' is negative. Although an equation similar to Eq. (11 1 was derived in Ref. 
[r3j for a general velocity-dependent force g{v\ this important issue was not reported and emphasis was only put 
on the additional entropy pumping contribution lm)t = —ASpu in the exponential factor (on the other hand, the 
time-reversal symmetry of the control force is discussed in the previous Ref. jl2)). It turns out however that changing 
the sign of 7' has also a significant consequence for the detailed FT, as discussed below. Hereafter, the stochastic 
process (dynamics) with 7' replaced by —7' is called the "backward" process (dynamics) for brevity [55]. 
Starting from Eq. ( 11 ), we now consider the ratio 



(12) 
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where the probability distributions po{xo,vo) and pi^xq^vq) for the initial and final states are still arbitrary at this 
stage. From Eq.(12), one readily obtains the integral fluctuation relation 



(13) 



where ^ ... ^ denotes a path integral average over all possible paths {xsj^gjo^tj from x{0) = xo,x{0) = vq to 
x{t) = xt,x{t) — vt, and is defined by 



<C w4[{a:s}] > = / dxodvo / dxtdvt 



(2:0, fo) 



VXsA[{Xs}]'P+[{Xs}\xo,Vo]po{xo,Vq) 



(14) 



for any trajectory-dependent functional ./[[{xs}]. As usual, one must make a suitable choice of the 'boundary' terms 
in Eq. ( 12 ) (those that only depend on the distributions of the initial and final states) to give a physical interpretation 
to the functional ^2?. The choice pi(a;(, uj) =p{x,v,t), where p{x , v , t) is the solution of the Kramers equation 

for the given initial distribution pq{xo,Vq), leads to 



R[{xs}] - AS^[{xs}] + AS+U = mxs}] , 



(15) 



which is the total entropy production in the super-system along the specific trajectory {xs}se[o,t]- Then Eq. (13) 
yields the integral fluctuation theorem (IFT) 



< e 



>= 1 



(16) 



already given in Ref. |13j . For a steady state at constant a characterized by the probability distribution 

Pss:a{x,v) = exp[-(f)aix,v)] , (17) 

the choice po(a;, w) =pi{x,v) = Pss;a{x,v) yields = Sss,q[{xs}], which from Eqs. Q-Q is given by 



T 



(18) 



where AE^ — (m/2)(W( — Vq) + [^^(xt) — VQ(a;o)] and A0q, = A5q, = 4'aixt,Vt) — (f'aixoTVo). In particular, the mean 
entropy production rate is 



<E,,,^[{a;,}] » ^ j' T-Teff 
t m T 



(19) 



where T^ff = m < (since T^ff < T for 7' > 0, the mean entropy production in the super-system is thus a 

positive quantity, as it must be). 

To derive the stronger detailed fluctuation theorem for the entropy production in a NESS, we now consider the 



probability that the functional S^^ q, [{x^}] given by Eq. ( 18 ) takes a specific value 1^ — at along the forward trajectory 

ftcr). This probability is defined by 



(to simplify the notation, the subscript a is dropped hereaf 

P+(I],,[{a:J]) = at) =« <5(I]«,[{a:,}] - at) » 



(20) 



Using Eq.(12), we then find 



dxodvo dxtdvt / 'Dxse^^^''''^^'P-[{xs}\xt,Vt]pss{xt,vt)S{T,ss[{xs}] - at) 



dxodvQ / dxfdvt 



'Dxs'P-[{xs}\xt,Vt]pss{xt,Vt)5{Y.ss[{xs}] - at) 



dxodvo dxtdvt / 'Dxs'P-[{xs}\xo,vo]pssixo,VQ)6(T.ss[{xs}] - at) (21) 

where the integration variables Xg and Xs, {xo,vo) and {xt,vt) have been interchanged to obtain the last equality. 
This equation can be written in the form of a detailed FT as 



P+{J:ss[{Xs}] ^ (Tt) 

P-{tss[{xs}] = ~<yt) 



(22) 
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where the functional ^^^[{xs}] is defined by 



Sss[{2;s}] = — Sss[{^s}] 



^ , ,/ 1 r , .2 t , 



(23) 



and P-CEssiixs}] = — is the probabihty that ^^^[{a;^}] takes the value — S along a trajectory in the backward 
process, given that the initial state is sampled from the steady-state probability Pssix,v) of the forward process. 

We thus see that the actual FT (which is valid for any length t of the trajectories) is more complicated that the 
one given in Ref. [T3]: the dynamics generating the stochastic trajectories in the numerator and the denominator are 
different, and so are the corresponding trajectory-dependent functionals (more precisely, the boundary terms —AE/T 
and in Eqs. (18 1 and (23) are identical whereas the sign of 7' is changed in the remaining 'bulk' term). Note also 



that Sss[{xs}] is not the entropy production functional in a steady state reached with the backward dynamics since 
the stationary distribution (and thus A^) is then different from the one given by Eq. (17) (as an example, see Eq. 
(26) below). In fact, as already pointed out, there is no stationary distribution with the backward dynamics if 7' is 
larger than the intrinsic friction 7 due to environment (7' > 7 is the current situation in a cold damping setup since 
the goal is to reduce the thermal noise as much as possible [71 [51 1^ ) . Nevertheless, the FT given by Eq. (22) is also 
valid in this case, as illustrated in Fig. 1 which shows the results of a numerical simulation of the Langevin equation 
for the harmonic potential studied in the next section. Note that the probability distributions for the forward and 
backward processes are quite different and that the latter (corresponding to a negative effective friction coefficient 
7 — 7') exhibits a long tail on the positive side. However, the relation ^^^(^^^[{a;^}] = at) ~ P-CSssHxs}] = —crt) e"'* 
is very well satisfied within the numerical accuracy of the calculation. 




FIG. 1: (Color online) Check of the detailed fluctuation theorem, Eq. ( |22[ ), for a stochastic harmonic oscillator with viscous 
dissipation and a velocity-dependent feedback force g{v) — —y'v. The figure shows the probability distribution functions 
P+(E3s [{a;^}] ~ crt) (green histogram) and P„(Esa[{a;3}] — at) (red histogram) for the forward and backward processes, 
respectively. The black dashed line represents the product P_(Ess[{a;s}] ~ —crt) e"^' obtained from the histogram. The model 
parameters in Eq. (24) are k — 0.2, 7 = 1, 7' = 1.2 (with m — 1 and T = 1), and the observation time is t = 2. The Langevin 
equation has been solved for 10'' realizations of the noise using Heun's method [27 with a time-step At — 0.001. 



It is clear that the complexity of the FT for a finite observation time t comes from the fact that the boundary 
and bulk terms in ^^^[{xs}] do not behave in the same way under time reversal and/or reversal of the feedback 
force. Therefore, one may expect some simplification if the contribution of the boundary term becomes negligible, 
which may occur in the long-time limit. This supposes, however, that 7' < 7 so that a steady state can be reached 
asymptotically in the backward process. Then Sss[{a;s}] becomes the actual entropy production in the steady state. 
This is illustrated in the next section by exact analytical calculations for the harmonic oscillator. 
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III. ILLUSTRATION ON THE STOCHASTIC HARMONIC OSCILLATOR 

A. Entropy production probability distribution in the steady state 

To illustrate the preceding discussion, we now consider the paradigmatic case of a harmonic potential and we 
calculate the probability distribution function (PDF) of the entropy production in a steady state. The stochastic 
harmonic oscillator with viscous dissipation is relevant to the dynamics of an AFM cantilever [5S], to the motion of a 
colloidal particle in an optical trap, and to many other practical applications or nano-electromechanical systems (see 
in particular Refs.[Tni [TT] for a recent application to the gravitational wave detector AURIGA). In general, it also 
permits a fully analytical analysis[5n]. In this case, Eq. ^ takes the very simple form 

mvs = ~kxs " (7 + l')vs + 6 (24) 
where k is the stiffness associated to the elastic force. The corresponding Kramers equation then writes 

dtp{x, V, t) = -vdxp{x, V, t) + —dy \ [kx + (7 + "i')v]p{x, v, t) + j—di,p{x, v,t)\ , (25) 

TO m 



which has the stationary solution 



^7 + 7' _r 7 + 7'r,, 2 , 21 



Pssix,v) = ^ exp{ — ^^[kx +mv ]} , (26) 



showing that the kinetic temperature of the Brownian system is T^ff — 7/(7 + 7')T. From Eq. (18), the total entropy 
production functional in the super-system is given by 

E,, [{x J] = Eii; (xo , «o , , «t ) + [{x J] (27) 

where 

r 

j:g\xo,vo,Xt,Vt) = ^[kix^t-xl)+m{v^-vl)] (28) 

and 

= fdsxl) (29) 

are the boundary and bulk contributions, respectively (note that both contributions vanish when 7' = since there 
is no other external force acting on the system which is then at equilibrium). I]ss[{xs}] is a quadratic functional of 
the noise and therefore its PDF is not Gaussian. Nevertheless, the generating or characteristic function defined by 

Z+(A,i)=<e-^^-[{"=>l > (30) 

can be explicitly calculated and the PDF is then recovered by taking the inverse Fourier transform 

P+(E,,[{xJ] =at) = — / dAZ+(A,t)e^'^* (31) 

where the integration is performed along the imaginary axis. We thus begin by computing Z-|_(A,t). 



Inserting Eq. (27) into Eq. (30), we obtain 



Z+{X,t) = / dxodvoPss{xo,vo) / dxtdvte-^^''<^^"'^o)+^^^i'-?,)] / Vx,r+[{xs}\xo,vo]e^ ^0 '^^ 

J J J(xq,Vo) 

(32) 

where 7^+[{j;s}|xo, t^o] is given by Eq. ([9]). This defines a new Lagrangian function 

Cx{xs,Xs,Xs) = -——(mxs + (7 + i)xs + kxsY + -Tr^xl , (33) 
47J 1 
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which can be rewritten as 
by introducing the A-dependent damping coefBcient 



1 2 T ^~ — T 
-^-;:^{mxs + ^{X)xs + kxs) — {k Xsis + mXsXs) 



7(A)^[(7 + 7')'-4A77T/' 



(34) 



(35) 



exp 



2m 



This yields 

r+[{x,}\xo.vo]e'^fod^-'i^) ^ 

where 

V^lixsjlxo^vo] = C exp 

Hence 



7 + 7' — 7 7 + 7' — 7 



7 



t- 



47T 



1 



2m 47r 



ds [mXs + 7X5 + kxg 



'P^[{xs}\xo,Vo] (36) 



(37) 



Z+{X,t) = e 



dxodvoPssixo,vo) / dxtdwtexp 



- ^ ^ ^Ij^ [H^t - ^0) + - «o)] 



X'[a::(s)] 'Pj[{xs}\xo,Vo] 
= e^^^^* / dxodvoPss{xo,vo) / dxtdvtcxp 



2mfi{X) + 4A7' 

4^r 



where yu(A) is defined by 



MA) . 7 + y-7(A)~2Ay 
2 m 



P^(a;t, ft, i|a;o, Wo, 0) 
(38) 



(39) 



and Pxf{xt,vt,t\xo,vo,0) is the transition probabihty (or propagator) corresponding to the damping coefficient 7(A). 
Note that Z^{X,t) as defined in the second line of Eq. (38) is properly normalized. Indeed, since 7(A = 0) = 7 + 7', 
one has /i(0) = and thus 



(40) 



Z+(0,i) EE< 1 > = / dxodvoPssixo,vo) J dxtdvtPj+Y{xt,Vt,t\xo,vo,0) 
dxodvo / dxtdvtP^+Y{xt,Vt,t;xo,vo,0) = 1 . 



Since the Langevin equation, Eq. (24), is linear and the noise is Gaussian, all stationary probability distributions 
are multivariate Gaussian distributions, and the explicit expression oi Pj{xt,Vt,t; xo,vo^O) is given by 



P=,{xt,vt,t; xo,vo,0) = ^^2(de!$,)i/2 ^^Pi-^^^^^ 
where $t is the matrix of time-correlation functions in the steady state 



(41) 



(bccAO) 0^.(0) 0..(t) (b,yit)\ /0XX(O) ^,^it) 0,,(t) \ 



_ I (l^vxiO) 4'vv{0) (t>vx{t) 4>vv{t) I _ 
4>xx{t) (j)vx{t) (l)xx{0) <l>xv{0) 
(l>xv{t) 4>vv{t) 0x1,(0) (j)vv{^) 



01,1,(0) (t)xx{-t) -(t^xxit) 

4>xx{t) (l)xx{~'t) (l)xx{'^) 
V kx[t) -^xx{t) 0,„(O) / 



and B is the 4-dimensional vector representing the initial and final conditions 



B 
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Here 4>xx{t) is the time-correlation function associated with an underdamped Langevin dynamics with damping 
coefBcient 7(A) [30], 



4>xx{t) 



-u:+\t\ 



-\t\ 



where 



7 ± y/y^~^-Akrr 



0J± = 



2m 



(42) 



(43) 



In particular, (pxxi^i) — lT/{'^k) and 4>vv{^) — lT/{'ym). 

Prom Eq. (41) we then compute the propagator _^(x(, wt, t|a;o, wo, 0) — P^jxt, Vt,t', xp, vp, Q)/p^{xq, vq) where 



p^{xp^vp) is obtained by replacing 7 + 7' by 7 in Eq. (26). Inserting into Eq. (38 1, we find 



where 



/ fc^(A) 



7T 





m^(A) 

k{^Ji{\) + 

, 2 x„/ 



V ^(^(A) + ^AY) / 

Carrying out the Gaussian integrals over xq, vq and Xt,Vt, we finally obtain the compact expression 

1 7+y *KA) 



Z+{X,t) 



[det(l + *tL)]i/2 -^(A) 



which is the main result of this section. 



(44) 



(45) 



Q 




FIG. 2: (Color online) Probability distribution P+(Esi,[{a;s}] = at) for k — 0.2, 7 = 1, 7' = 0.2 (green histogram) and 7' = 1.2 
(red histogram). The observation time is t = 2. The numerical inverse Fourier transforms of Eq. (45 I (shown as dashed black 



lines) are compared to the histograms of Sss[{a;s}] obtained from the simulation of the Langevin equation. 



This generating function is a complicated function of A and its inverse Fourier transform can only be performed 
numericall y. N ote that it was implicitly assumed in the above calculations that the damping coefficient 7(A) is real. 
From Eq. (35), this implies that A < Xmax = (7 + 7')^/(477') on the real axis. On the other hand, the integration 
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in Eq. (31) is performed along the imaginary axis and therefore quantities hke P;y[{a;s}|a:o, I'd] (defined by Eq. (37)) 



become complex and loose their physical meaning. This may cast some doubt on the validity of the final result, Eq 
(45). However, one can show that the calculation remains valid and that the inverse Fourier transform of Eq, 

}] ~ at) (this will be c 



(45) 



ear 



is indeed a real quantity which correctly defines the probability distribution P_)-(E[{a: 

below when considering the asymptotic long-time behavior). As shown in Fig. [2j the numerical Fourier transform of 
the theoretical expression is indeed in excellent agreement with the histograms of S]+[{a;s}] obtained from the direct 
simulation of the stochastic process, for both 7' < 7 and 7' > 7. One can also clearly observe that the distributions 
are non-Gaussian. 

Since 7(0) =7 + 7' and ^(0) = (which imphes L = 0), one readily sees from Eq. (45) that ^+(0, t) is properly 
normalized. On the other hand, it is not immediately apparent that the integral fluctuation theorem Z^{l,t) = 1 is 
satisfied. One needs to distinguish the two cases 7 > 7' and 7 < 7'. In the first case, one has 7(1) =7 — 7' so that 
/i(l) = 0. The calculation of the determinant in Eq. (45 1 then gives 



Vdet(l + *tL) = 



1 + 1' 
7-7' 



(46) 



and thus Z+{l,t) = 1, as it must be. Note that this result is obtained without using the explicit expressions of 
the time-correlation functions, but only their values at t — 0. In the second case, one has 7(1) = 7' — 7 and 
/i(l) = (7 — 7')/m, and the calculation of the determinant gives 

v/dct(l + *tL) = ^""^^^^'T^ \ cb^At)q^,,{t) - </),„(O0.x(O] 
7^J ^ 



(w-i- + a;_)•^w_|_a;_ 



7 + 7 i-Y f 
-e 



7' - 7 



(47) 



Inserting into Eq. (45) yields the correct result Z+{l,t) = 1. Remarkably, in this case, we had to use the explicit 
expressions of the time-correlation functions. 



In order to check the detailed FT expressed by Eq. ( 22 ) , one needs to calculate the generating function Z_ (A, t) of the 
functional Ess[{xs}] in the backward process. Formally, one can follow the same steps as in the preceding calculation, 
at least up to Eq. (38 ) (replacing 7' by —7'). To proceed further, however, one needs to compute Pxf{xt,Vt, t; xq, wo, 0), 
that is to solve the Kramers equation for the backward process with ^^^^(aJo, 'I'o) as initial condition (that is with the 
stationary distribution of the forward process) . This is a complicated calculation which we have not performed (see 
Fig. 1 for a numerical check of Eq. (22 )), and in the following we shall only consider the asymptotic long-time regime. 



We just note that if the conventional FT were to hold exactly one would have 



(48) 



From Eq. ([35|, one has 7+(l — A) = 7-(A) so that /i+(l — A) = /i-(A) where the indices -I- and — refer to 7' and 
—7', respectively. Therefore one also has #t^+(l — A) = $t._(A), assuming that the backward process has reached a 
steady-state (which implies that 7' < 7). The only function of A that does not have a simple symmetry is the matrix 
L, which is not surprising since it contains the information about the initial conditions. 



B. Long-time behavior and large deviation function 

To complete this study we now consider the behavior of P+(Sss[{a;s}] = crt) for t much larger than the effective 
viscous relaxation time = m/(7 -f 7'). We expect a large deviation form |31j 

P+(S,,[K}] =at)^e''(^)* (49) 
where h{a) is the large deviation function (LDF) defined by 

h{a) = lim \ lnP+(E,,[{xJ] - nt) . (50) 



As usual, the large-^ behavior of the PDF can be extracted from the integral representation (31 ) by using a saddle- 
point approximation and taking care of the possible presence of singularities in the integrand. Hence we first need 
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to derive the asymptotic form of Z+{X,t) from Eq. 



(45), which is easily done by observing that the real parts 
= (7 + 7')^/(477'). The time-correlation functions 



of (jj± are always positive if 7(A) is real, that is if A < Xmax 
4'xx{i)i 4'xv{i)^ 4>vv{i) then go to zero as t — )> cxd and the matrix $j becomes diagonal. The determinant of 1 + in 
Eq. (45) is readily calculated, and we obtain 



Z+(A,t)^g(A)e*^(^) 



with 



5(A) 



4(7 + i)i 

|(7 + 7'+7)2"4A27'2| 



(51) 



(52) 



and /i(A) defined by Eq. (39). (Note that this calculation is not valid for A = 1 in the case 7' > 7 since then 
\/det(l + $tL) — 0, as can be seen from Eq. (47). In this case Z+{l,t) = 1 for all t.) 

By definition, ^(A) = hm(_i.oo \nZ{X,t)/t is the cumulant generating function [3T] and the saddle point A*(cr) is then 
solution of the equation 



with 



m'(a) 







7' 7 -7(A) 



m 



7(A) 



(53) 



(54) 



from Eq. (39). Since lim 



.m'(a) 



On the other hand, for a < j' /m, the solution of Eq. ( |53[ ) is given by 

[7'^ - ma{j + 7')] [7'^ - TOcr(7 + 7') 



A* (a) 



"f' /m, we see that the saddle point equation has no solution for a > /m 

277'] 



477' (7' 



(55) 



which is a function of a that monotonically decreases from Xmax to —00 as a increases from —00 to j'/m [32j . 
However, we also note from Eq. (52) that the prefactor g{X) diverges when 27'A = ±(7 + 7' + 7), and we thus have 
to determine for which values of a the saddle point coalesces with a pole of g{X). To proceed further, we need to 
consider the two cases 7 > 7' and 7 < 7' separately. 



1. 7 > 7 

In this case g{X) has a simple pole located on the real axis at A = X„iin = ~(27-l-7')/7' ^-nd by solving the equation 
A* (a) — Xmin we find that the saddle point hits this pole at cr = ct* with 

7' 27 + 7' 

-7 = — o , , • (56) 
m 37 + 7 

Starting from a ~ —00, one can thus safely deform the contour of integration through the saddle point as long as 
a < a* . The LDF is then given by the Legendre transform of the cumulant generating function|31j 

h{(j)= pl{X*) + X*(j (57) 

which yields 

M.).>..w=- i^;-7'7-^^'f . (58) 

4m77'(7' — ma) 

On the other hand, for a > a* , the steepest-descent contour must cross the pole and the leading contribution to the 
integral comes from the pole (see the appendix for more details). The LDF is then a linear function of a 

h{a) = n{Xmin) + XrmnCr (59) 

which yields 

h[a) = h2((7) = , — cr . (60) 

m 7' 
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FIG. 3: (Color online) Large deviation function h{a) for 7 > 7' (m = 1,7 = 1,7' = 0.2). The vertical dashed line marks 
the position of a* that separates the branches hi (a) and /i2 (o") . The (red) dotted-dashed line represents the branch hi {a) for 
a > a* which diverges a,t a — 7'/m. Note that hi{a) — for a — j'^ /[m{'Y + 7')], the mean entropy production rate. 



The behavior of the LDF as a function of cr is illustrated numerically in Fig. [s] (one can easily check that h{a) 
and its first derivative are continuous at ct = ct*). In the appendix, we give the complete asymptotic form of the 
probability distribution P_|_(I]ss[{a;s}] = at), taking into account the proximity of the saddle point to the pole as 
explained in Ref. [33], and the contribution of the residue of g{X) at A = A„ii„ when the pole is crossed. Interestingly, 
the PDF in the long-time limit becomes independent of k, the stiffness of the harmonic oscillator (see Ref. |M| for a 
similar observation). 

Before considering the case 7 < 7', let us briefly comment these results. We first note that the pole in g{X), 
which limits the position of the saddle point, appears when the average over the initial and final states is perfor med . 



Its presence can be traced back to the contribution of the boundary term in Sss[{a;s}] given by Eq. (28 1 



The singularity arises because the position and the velocity of the particle are unbounded and Gaussian distributed 



according to the steady-state probability distribution ( |26[ ). As a consequence, the PDF of Y^i]) has an exponentially 
decreasing tail (specifically, of the form e ^' ''^'') and large fluctuations of order t may occur which cannot be 



neglected in the sum + Y,si' despite the fact that Ess'' is not extensive in time. On the other hand, there is no 
singularity if the initial and final positions and velocities are fixed. Such an interplay between boundary and bulk 
terms is well documented in the literature on large deviationsjM] and fluctuation relations (see e.g. Refs. O I35H40] 
and more recently Ref. [34j ) . In the present case, however, this interplay is a bit unusual. Indeed, the slope of h2{cr) 
is not equal to — (7 + 7')/7'; which means that it is not simply imposed by the exponential tail of the PDF of Els'* (in 
contrast, for instance, with the model studied in Ref. 0^). Clearly, one cannot treat the fluctuations of the boundary 
and bulk contributions independently, even asymptotically. This is all the more remarkable that the latter (divided 
by t) is bounded by /m as can be readily seen from Eq. (29). This is actually the origin of the divergence in A*(ct) 
and hi{a) at cr = 7'/m. However, since a* < /m, this divergence occurs in the region where the large fluctuations 
are described by ft.2(o') and it is thus harmless (see Fig. |3|. 



.(2) 



^(1) 



We can now come back to the detailed fluctuation theorem (Eq. (22)) by noting that the function hi{a) possesses 
the symmetry 



h^{cr) — (—0") = cr 



(61) 



where the superscripts -I- and 
asymptotic limit of the FT 



refer to 7' and —7' respectively (hence h^{cr) = hi{a)). This is precisely the 



^. 1 P+{J:ss[{Xs} = <7t] 

lim - In ^ = cr 

t^^t P_(I],,[{xJ = ~at] 



(62) 
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provided that the large deviation forms 



(63a) 
(63b) 



are both vahd. We already know that Eq. (63a) is valid for a < a* . Eq. (63b) is also valid when the fluctuations 
of the boundary term becomes negligible so that it is irrelevant to sample the initial state of the backward process 
with the steady-state probability of the forward process. Then Eq. ( |5l| ) (with 7' replaced by —7' ) correctly describes 
asymptotically the generating function Z-{X,t) of the functional Sss[{2;s}]. This occurs for a > fil = fT*(— 7'). One 
can check from Eq. (561 that a* > —a*_ and therefore Eq. (61 1 is indeed the asymptotic expression of the detailed 
FT for o- < -crl (in tEe case displayed in Fig. [3| cr* = 0.1375 and a*_ « -0.129). 



2. 7 < 7' 

We now turn to the physically more relevant case 7 < 7'. The new feature is the presence of a second pole in g{X) 



at A = 1 (see also the remark after Eq. (52)). By solving the equation A*((t) = 1, we find that the saddle point and 
the pole coalesce at a — a** with 

a- ^ . (64) 

m 7—7 

This value is smaller than cr* and therefore the LDF is described by the function /ii(cr) (the Legendre transform of the 
cumulant generating function) in the interval [(T**,(T*] only. On the other hand, for a < a**, the leading contribution 
to the integral comes from the pole at A = 1 and the LDF is again linear 

h{(7) = fi{l) + a (65) 

which yields 

h{a) = hsia) ^ -^'^^ + a . (66) 
m 

Finally, for cr > cr*, the contribution from the other pole is dominant and h(a) = /12 (cr) like in the case 7 > 7'. 




FIG. 4: (Color online) Large deviation function h{a) for 7 < 7' (m = 1,7 = 1,7' = 5). The vertical dashed lines mark the 
positions of cr** and cr* that separate the branches ft3(cr), /ii (cr), and ^2(0"). The (red) dotted-dashed line represents the branch 
hi (a) for a > a* which diverges at cr = 7'/m. 
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The behavior of the LDF as a function of cr is illustrated numerically in Fig. [4] and the complete asymptotic 
expression of the PDF (taking into account the presence of the two poles in g{X)) is also given in the appendix. Fig. 
[5] in the appendix confirms that this expression correctly describes the PDF when the observation time t is very large. 

Note that a* ~ cr** ~ j' /m when 7'/7 ^ 1 so that the domain of validity of the central branch ft.i(cr) becomes very 
small. In any case, the symmetry relation (61 1 cannot be interpreted as the asymptotic expression of the detailed FT 
when 'y' > 7. Indeed, as already stressed, the system does not reach a steady state with the 'backward' dynamics and 
Eq. (45 1 does not describe the generating function Z_{X, t) of the functional ^^^[{xs}] (we recall that the expressions 
of the steady-state time-correlation functions have been used to derive this equation). 

Finally, let us stress that a similar analysis can be performed for the probability distribution functions of the heat 
adsorbed by the oscillator or the injected power (the expressions of the characteristic functions are quite similar to 
Eq. (45)). The analytical results can be directly compared to the data collected by the gravitational wave detector 
AURIGA which are given in Ref.|llj. 



IV. CONCLUDING REMARKS 



To summarize, we have revisited the model of a classical molecular refrigerator described by an underdamped 
Langevin equation with a feedback force proportional to the velocity. Unlike the viscous force due to the environment, 
the feedback can be seen as a virtual viscous force that creates dissipation without introducing fluctuations. This 
modifies the entropy production in the system and the contribution of the feedback mechanism to the entropy must 
be included in the second law and the fluctuation theorems, as discussed previously 1121 . However, we have shown 
that the detailed fluctuation theorem has a more complicated interpretation than originally suggested [T3]. This results 
from the fact that the sign of 7', the friction coefRcient associated to the feedback force, must be changed in order to 
determine the appropriate backward (time-reversed) path corresponding to a given forward path in the path integral 
approach. This kind of issue has already been discussed in the literature j24j but it takes a special importance here due 
to the friction-like character of the feedback force[l2]. For instance, this implies that the system is heated and cannot 
reach a stationary state in the backward process if 7' is larger than 7, the intrinsic friction due to the environment. 
7' > 7 is in fact the common experimental situation. By solving analytically the harmonic oscillator and computing 
the probability distribution of the total entropy production in a NESS we have shown that the regime of fluctuations 
in the cooling mode (the usual forward process) also depends on whether the ratio 7'/7 is smaller or larger than 1. In 
particular, the large time behavior of the PDF, as described by the large deviation function, is controlled by a subtle 
and rather unusual interplay between boundary and bulk contributions. This is a remarkable feature taking into 
account the simplicity of the model and it might be an interesting challenge to check this behavior experimentally. 
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Appendix A: Asymptotic expression of the PDF 

In this appendix we give the asymptotic form of the probability distribution P+(I]ss[{xs}] = crt) in the long-time 
limit. To take into account the proximity of the saddle point to a pole in g{X) (for instance, for 7 > 7', the pole at 
Xmin which is reached when ct = cr*), we write 

g(A)= ,^-; +.g(A) (Al) 

where glj^ is the residue of g{X) at Xmin and g{X) is the regular part. We then treat the two contributions to the 
contour integral as explained in Ref.[33| and add the contribution of the residue when the contour has to cross the 
pole (see also Ref . |34) for a similar calculation). Similarly, in presence of the other pole at A = 1 for 7 < 7', we write 

g(A)^ , +^+9W ■ (A2) 

A — Amin A i 



Skipping the details, we find: 
a) For 7 > 7', 
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p+{j:ss[{xs}]^<jt) « <^ 



g/ll(<T)t 



for a < a* 



/2(a) + e''^(-)*5li[l - ierfc(v/M^)] for <7* < a < ^ 



(A3) 



* „/l2(o-)t 



for CT > ^ 



where 



u{(7) = ft,2(o-) - ft-l(o-) , 



(A4) 



g(A*(a)) 9-1 

g(A*(a)) , 9U 
v/2/i"(A*(a)) 2vAIM ' 



(A5) 



and 



5-1 



(7 + 7')(37 + 7')' 



27'(27 + 7') 



'^2 



(A6) 



b) For 7 < 7', 




FIG. 5: (Color online) PDF of the entropy production for 7 = 1 and 7' = 5 (T = 1, m = 1, it = 0.2) and t = 20, 30, 40 (from 
top to bottom). The solid black lines represent the analytical asymptotic expressions given by Eqs. (A7l and the dashed blue 
lines are the corresponding numerical inverse Fourier transforms of Z+{X,t). The red points are obtained from the numerical 
simulation of the Langevin equation for t — 20. The vertical dashed lines mark the positions of a* (right) and a** (left). 
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( f,hi{a-)t 



P+{E,s[{Xs}]=<jt) « < 



TTt 



'lit 
g*_^eh2{a)t 



where 



for a < a** 
for (7** < cr < cr* 
for cr* < CT < ^ 

— — m 

for a > ^ 
(A7) 



v{<7) = h3{a) - hi{a) 



(A8) 



and 



.9(A*(<t)) 


.9-1 


.9-1 


v/2M"(A*(a)) 


2v/M(fT) 


2V«M 


5(A*(a)) 




^ .9-1 


2v'2^"(A*(a)) 




2V^(a) 


5(A*(a)) 




^ .9" 


2V2A*"(A*(a)) 







5-1 = 



(Y-7)^(Y + 7) 
27'3 



(A9) 



(AlO) 



As shown in Fig. [5j the above asymptotic expressions are in excellent agreement with the numerical inverse Fourier 
transform of Z+(A,f) (Eq. (45 1). In particular, we note that the small discrepancies on the right hand side (for 
a > a*) diminish as t increases. For t = 20, there is also a good agreement with the numerical simulation of the 
Langevin equation. 
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